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Abstract 

We present a generalized version of the ITIM algorithm for the identification of interfacial molecules, 
which is able to treat arbitrarily shaped interfaces. The algorithm exploits the similarities between 
the concept of probe sphere used in ITIM and the circumspherc criterion used in the a-shapes 
approach, and can be regarded either as a reference-frame independent version of the former, or 
as an extended version of the latter that includes the atomic excluded volume. The new algorithm 
is applied to compute the intrinsic orientational order parameters of water around a DPC and a 
cholic acid micelle in aqueous environment, and to the identification of solvent-reachable sites in 
four model structTires for soot. The additional algorithm introduced for the calculation of intrinsic 
density profiles in arbitrary geometries proved to be extremely useful also for planar interfaces, as 
it allows to solve the paradox of smeared intrinsic profiles far from the interface. 
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I. INTRODUCTION 



Capillary waves represent a conceptual problem for the interpretation of the properties of 
liquid-liquid or liquid-vapor planar interfaces, because long-wave fluctuations are smearing 
the density profile across the interface and all other quantities associated to it. This is 
usually overcome by calculating the density profile using a local, instantaneous reference 
frame located at the interface, commonly referred to as the intrinsic density profile, p{z) = 
{A~^ S {z — Zi + ^{xi, Hi))), where {xi,yi,Zi) is the position of the i-th atom or molecule, 
and the local elevation of the surface is assuming the macroscopic surface normal 

being aligned with the Z axis of a simulation box with cross section area A. During the 
last decade several numerical methods have been proposed to compute the intrinsic density 
profiles at interfaces'"^. Despite several differences in these approaches, they are, in general, 
providing consistent distributions of interfacial atoms or molecules^ and density profiles'^. 
Among these methods, ITIM^ proved to be an excellent compromise between computational 
cost and accuracy^, but it is limited to macroscopically flat interfaces, therefore there is a 
need to generalize it to arbitrary interfacial shapes. 

Before these works, albeit for other purposes, several surface-recognition algorithms have 
been devised, and will be briefly mentioned below. All of them are possible starting points 
for the sought generalization under the condition that, once applied to the special case of a 
planar interface, they lead to consistent results with existing algorithms for the determination 
of intrinsic profiles. 

Historically, the first class of algorithms addressing the problem of identifying surfaces 
was developed to determine molecular areas and volumes. The study of solvation proper- 
ties of molecules and macromolecules (usually, proteins) might require the identification of 
molecular pockets, or the calculation of the solvent-accessible surface area for implicit solva- 
tion models^. Two intuitive concepts are commonly used to describe the surface properties 
of molecules, namely, that of solvent-accessible surface^''° (SAS), and that of molecular 
surface''"'^ (MS, also known as solvent excluded surface, or Connolly surface). The MS 
can be thought as the surface obtained by letting a hard sphere roll at close contact with 
the atoms of the molecule, to generate a smooth surface made of a connection of pieces of 
spheres and tori, which represents the part of the van der Waals surface exposed to the 
solvent. During the process of determining the surface, interfacial atoms can be identified 
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using a simple geometrical criterion. Many approximated^^"^ or analyticaP^'^^'^^"^^ meth- 
ods have been developed to compute the MS or the SAS. In general, these methods are 
based on discretization or tessellation procedures, requiring therefore the determination of 
the geometrical structure of the molecule. Other methods which allow to identify molecular 
surfaces include the approaches of Willard and Chandler^ or the Circular Variance method 
of Mezei^^ . Incidentally, the way the MS is computed in the early work of Greer and Bush^^ 
resembles very closely the ITIM algorithm^. 

Prom the late 1970s, the problem of shape identification had started being addressed 
by a newly born discipline, computational geometry. In this different framework, several 
algorithms have been actively pursued to provide a workable definition of surface, and in 
particular the concept of a-shapes^^'^^ showed direct implications for the determination of 
the molecular surfaces^'^'^^. The approach based on a-shapes is particularly appealing due to 
its generality and ability to describe, besides the geometry, also the intermolecular topology 
of the system. 

Noticeably, none of these methods - to the best of our knowledge - has ever been em- 
ployed for the determination of intrinsic properties at liquid-liquid or liquid-gas interfaces. 
Prompted by the apparent similarities between the usage of the circumspherc in the alpha 
shapes and that of the probe sphere in the ITIM method, as we will describe in the next 
section, we investigated in more detail the connection between these two algorithms. As a 
result, we developed a generalized version of ITIM (gitim) based on the a-shapes algorithm. 
The new GiTiM method consistently reproduces the results of ITIM in the planar case while 
retaining the ability to describe arbitrarily shaped surfaces. In the following we describe 
briefiy the alpha shapes and the ITIM algorithms, explain in detail the generalization of the 
latter to arbitrarily shaped surfaces, and present several applications. 

II. ALPHA SHAPES AND THE GENERALIZED ITIM ALGORITHM 

The concept of a-shapes was introduced several decades ago by Edelsbrunner^^'^^. To 
date the method is applied in computer graphics application for digital shape sampling 

and processing, in pattern recognition algorithms and in structural molecular biology^®. The 
starting point in the determination of the surface of a set of points in the a-shapes algorithm 
is the calculation of the Delaunay triangulation, one of the most fruitful concepts for compu- 
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tational geometry^ which can be defined in several equivalent ways, for example, as the 
triangulation that maximizes the smallest angle of all triangles, or the triangulation of the 
centers of neighboring Voronoi cells. The idea behind the a-shapes algorithm is to perform 
a Delaunay triangulation of a set of points, and then generate the so-called a-complex from 
the union of all k-simplices (segments, triangles and tetrahedra, for the simplex dimension 
k—1,2 and 3, respectively), characterized by a k-circumsphere radius (which is the length 
of the segment, the radius of the circumcircle and the radius of the circumsphere for k—1,2 
and 3, respectively) smaller than a given value, a (hence the name). The a-shape is then 
defined as the border of the a-complex, and is a polytope which can be, in general, concave, 
topologically disconnected, and composed of patches of triangles, strings of edges and even 
sets of isolated points. In a pictorial way, one can imagine the a-shape procedure as growing 
probe spheres at every point in space until they touch the nearest four atoms. These spheres 
will have, in general, different radii. Those atoms that are touched by spheres with radii 
larger than the predefined value a are considered to be at the surface. 

An example of the result of the a-shapes algorithm in two dimensions is sketched in 
Fig. la. The ITIM algorithm is based instead on the idea of selecting those atoms of one 
phase that can be reached by a probe sphere with fixed radius streaming from the other 
phase along a straight line, perpendicular to the macroscopic surface. An atom is considered 
to be reached by the probe sphere if the two can come at a distance equal to the sum of 
the probe sphere and Lennard- Jones radii, and no other atom was touched before along the 
trajectory of the probe sphere. In practice, one selects a finite number of streamlines, and 
if the space between them is considerably smaller than the typical Lennard- Jones radius 
Rp, the result of the algorithm is practically independent of the location and density of the 
streamlines. The same is not true regarding the orientation of the streamlines; this is a direct 
consequence of the algorithm being designed for planar surfaces only. The basic idea behind 
the ITIM algorithm are sketched in Fig. lb. A closer inspection reveals that the condition of 
being a surface atom for the ITIM algorithm resembles very much that of the a-shapes case. 
Quadruplets of surface atoms identified by the ITIM algorithm have the characteristic of 
sharing a common touching sphere having the same radius as the probe sphere. In this way, 
one can see the analogy with the a-shapes algorithm, the Rp parameter being used instead 
of a. The most important differences in the a-shapes algorithm with respect to ITIM are the 
absence of a volume associated with the atoms, and its independence from any reference 
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frame. We devised, therefore, a variant of the a-shapes algorithm that takes into account 
the excluded volume of the atoms. 

In the approach presented here the usual Delaunay triangulation is performed, but the 
a-complex is computed substituting the concept of the circumsphere radius with that of the 
radius of the touching sphere, thus introducing the excluded volume in the calculation of 
the a-complex. Note that this is different from other approaches that are trying to mimic 
the presence of excluded volume at a more fundamental level, like the weighted a-shapes 
algorithm, which uses the so-called regular triangulation instead of the Delaunay one^^. In 
addition, in order to eliminate all those complexes, such as strings of segments or isolated 
points, which are rightful elements of the shape, but do not allow a satisfactory definition 
of a surface, the search for elements of the a-complex stops in our algorithm at the level 
of tetrahedra, and triangles and segments are not checked. In this sense GITIM can provide 
substantially different results from the original a-shapes algorithm. 

The equivalent of the ct-complex is then realized by selecting the tetrahedra from the 
Delaunay triangulation whose touching sphere is smaller than a probe sphere of radius Rp, 
and the equivalent of the a-shape is just its border, as in the original a-shapes algorithm. 
The procedure to compute the touching sphere radius is described in the Appendix. 

In the implementation presented here, in order to compute efficiently the Delaunay tri- 
angulation, we have made use of the quickhuU algorithm, which takes advantage of the 
fact that a Delaunay triangulation in d dimensions can be obtained from the ridges of the 
lower convex hull in d+ 1 dimensions of the same set of points lifted to a paraboloid in the 
ancillary dimension^^. The quickhuU algorithm employed here^° has the particularly advan- 
tageous scaling 0(iVlog(z/)) of its computing time with the number N and u of input points 
and output vertices, respectively. 

A separate issue is represented by the calculation of the intrinsic profiles (whether profiles 
of mass density or of any other quantity) as the distance of an atom in the phase of interest 
from the surface is not calculated as straightforwardly as in the respective non-intrinsic 
versions. For each atom in the phase, in fact, three atoms among the interfacial ones have 
to be identified in order to determine by triangulation^ the instantaneous, local position 
of the interface. This issue will be discussed in Sec. Ill for the planar, for the spherical or 
quasi-spherical and for the general case: here we simply note that we turned down an early 
implementation of the algorithm that searches for these surface atoms, based on the sorting 
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of the distances using 0(A'"log A^) algorithms hke quicksort, in favor of a better performing 
approach, based on kd-trees^^'^^, a generahzation of the one-dimensional binary tree, which 
are still built in a C(A'"logA'") time, but allow for range search in (typically) C(logiV) time. 

III. COMPARISON BETWEEN THE ITIM AND THE GITIM METHODS 

We have compared the results of the ITIM and GiTiM algorithms applied to the wa- 
ter/carbon tetrachloride interface composed of 6626 water and 966 CCI4 molecules. The 
water and CCI4 molecules have been described by the TIP4P model^^, and by the potential 
of McDonald and coworkers^^, respectively. The molecules have been kept rigid using the 
SHAKE algorithm^^. This simulation, as well as the others reported in this work have been 
performed using the Gromacs^® simulation package employing an integration time step of 
1 fs, periodic boundary conditions, a cutoff at 0.8 nm for Lennard- Jones interactions and 
the smooth Particle Mesh Ewald algorithm^^ for computing the electrostatic interaction, 
with a mesh spacing of 0.12 nm (also with a cut-off at 0.8 nm for the real-space part of the 
interaction) . All simulations were performed in the canonical ensemble at a temperature of 
300K using the Nose-Hoover thermostat ^^'^^ with a relaxation time of 0.1 ps. A simulation 
snapshot of the H2O/CCI4 interface is presented in Fig. 2, where the surface atoms identi- 
fied by the gitim algorithm using a probe sphere radius of 0.25 nm are highlighted using a 
spherical halo. 

We have used the itim and GiTiM algorithms to identify the interfacial atoms of the 
water phase in the system, for different sizes of the probe sphere. In general, GiTiM identifies 
systematically a larger number of interfacial atoms than ITIM for the same value of the probe 
sphere radius as it is clearly seen in Fig. 3. Remarkably, for values of the probe sphere 
radius smaller than about 0.2 nm (compare, for example, with the optimal ITIM parameter 
Rp — 0.125 nm suggested in Ref. 6), the interfacial atoms identified by GiTiM show the onset 
of percolation. The reason for this behavior traces back to the fact that ITIM is unable to 
identify voids buried in the middle of the phase, as it is effectively probing only the cross 
section of the voids along the direction of the streamlines. This difference could explain the 
higher number of surface atoms identified by GITIM, as voids in a region with high local 
curvature (or, in other words, with a local surface normal which deviates significantly from 
the macroscopic one) will not be identified as such by ITIM. In GiTiM, on the contrary, probe 
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spheres can be thought as inflating at every point in space instead of moving down the 
streamhnes, and this is the reason why the algorithm is able to identify also small pockets 
inside the opposite phase. 

It is possible to make a rough but enlightening analytical estimate of the probability for 
a probe sphere of null radius in the ITIM algorithm to penetrate for a distance C in a fluid of 
hard spheres with diameter a and number density p. Using the very crude approximation of 
randomly distributed spheres, the probability po to pass the first molecular layer, at a depth 
C = cr is the effective cross section ■pQ — \ — '^p^l'^a^^ and that of reaching a generic depth 
C can be approximated as = Po*^"^' where k = ln(l/j9o)/c defines a penetration depth. 
Therefore, using a probe sphere with a null radius, ITIM will identify a (diffuse) surface at 
a depth while gitim will identify every atom as a surface one. For water at ambient 
conditions, the penetration is ~ 0.186 nm, a distance smaller than the size of a water 
molecule itself. This could explain why in Ref. 6, even using a probe sphere radius as small as 
0.05 nm, almost only water molecules in the first layer were identified as interfacial ones by 
ITIM (see the almost perfectly Gaussian distribution of interfacial water molecules in Fig.9 
of Ref. 6). 

Nevertheless, it is important for practical reasons to be able to match the outcome of 
both algorithms. It turns out that choosing Rp so that the average number of interfacial 
atoms identified by both algorithms is roughly the same leads also, not surprisingly, to very 
similar distributions. The probe sphere radius required for GiTiM to obtain a similar average 
number of surface atom as in ITIM can be obtained by an interpolation of the values reported 
in Fig. 3. An example showing explicitly the interfacial atoms identified by the two methods 
{Rp = 0.2 nm for ITIM and Rp = 0.25 nm for GITIM) is presented in Fig. 4: roughly 85% 
of surface atoms are identified simultaneously by both methods, demonstrating the good 
agreement between the two methods once the probe sphere radius has been re-gauged. The 
condition of identifying the same atoms as interfacial ones is much more strict that any 
condition on average quantities, like the spatial distribution of interfacial atoms or intrinsic 
density profiles. Hence, it is expected that a good agreement on such quantities can also be 
achieved. 

The intrinsic density profiles of water and carbon tetrachloride are reported in Fig. 5, as 
computed by ITIM and GiTiM, respectively, with the interfacial water molecules as reference. 
The procedure for identifying the local distance of an atom from the surface is in its essence 
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the same as described in Ref. 7. Starting from the projection Pq — {x, y) of the position of 
the given atom to the macroscopic interface plane, the two interfacial atoms closest to Pq 
are found (their position on the interface plane being Pi and P2, respectively). The third 
closest atom with projection P3 has then to be found, with the condition that the triangle 
P1P2-P3 contains the point Pq- A linear interpolation of the elevation of Pq from those of 
the other points is eventually performed, and employed to compute the distance z — ^{x, y) 
which is used to compute the intrinsic density profile. 

Efficient neighbor search for the Pi, P2 and candidate P3 atoms is implemented using 
kd-trees^^ as discussed before. The two pairs of profiles are very similar, besides a small 
difference in the position and height of the main peak of the CCI4 profile (curves on the 
right in Fig. 5) and in the minimum of the water profile (curves on the left in Fig. 5) right next 
to the surface position, which are anyway compatible with the differences observed between 
various methods for the calculation of intrinsic density profiles^. The delta- like contribution 
of the water molecules at the surface is included in the plot in Fig. 5, and defines the origin 
of the reference system. Negative values of the signed distance from the interface correspond 
to the aqueous phase. 

IV. THE PROBLEM OF NORMALIZATION OF DENSITY PROFILES 

Before applying GITIM to non-planar interfaces, one important issue has still to be solved, 
namely that of the proper calculation of intrinsic density profiles in non- planar geometries. 
In general, one uses one-dimensional density profiles (intrinsic or non-intrinsic) when the 
system is, or is assumed to be, invariant under displacements along the interface, so that the 
orthogonal degrees of freedom can be integrated out. When the interface has a non-planar 
shape, one needs to use a different coordinate system. In the case of a quasi-spherical object 
for example, one could use the spherical coordinate system to compute the non-intrinsic 
density profile, and normalize each bin by the integral of the Jacobian determinant, that is 
the volume of the shell at constant distance from the origin. In the intrinsic case, however, 
it is necessary to know at every time step the volume of the shells at constant distance from 
the interface. 

The volume of shells at constant intrinsic distance can, in principle, be calculated at each 
frame by regular numerical integration, but this would require a large computing time and 
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storage overhead. Here, instead, we propose to employ an approach based on simple Monte 
Carlo integration: in parallel with the calculation of the histograms for the various phases, 
we compute also that of a random distribution of points, equal in number to the total atoms 
in the simulation. The volume of a shell can be estimated as box volume times the ratio 
of the number of points found at a given distance and the total number of random points 
drawn. We are following the heuristic idea that for each frame j one does not need to know 
the volume of the shell Vj{r) with a precision higher than that of the average number of 
atoms in it, Nj{r). In addition, we assume that the surface area of the interface is large 
enough for the shell volume variations SVj{r) to be small with respect to its average value 
V{r) = Vj{r)/N. The average density 



When the relative volume changes |(5V^/V^| are small, one can therefore simply normalize the 



The correctness of our assumption is demonstrated incidentally by the application of 
this normalization once again to the planar case. The thin lines in Fig. 5 represent the 
ITIM intrinsic mass density profile of water and carbon tetrachloride, using the Monte Carlo 
normalization scheme instead of the usual normalization with box cross sectional area and 
slab width. Close to the interface, the Monte Carlo normahzation gives results which are 
fully compatible with the usual method, showing that the accuracy of the volume estimate is 
adequate. On the other hand one can see that far from the interface the two profiles behave 
quite differently. The case with usual iiormalizatioii decays slowly to zero: this effect is due to 
the presence of the second interface, whose profile is smeared again by capillary waves. The 
case with Monte Carlo normalization, on the contrary, shows that it is possible to recover 
the proper intrinsic density also at larger distances, and features such as the fourth peak 
at 2 nm, which are completely hidden in the normal picture, can be revealed. This shows 
that the use of the proper, curvilinear coordinate system is of fundamental importance also 




(1) 



can be approximated as 





histogram N{r) — Nj{r)/N by the average volume V{r) obtained by the Monte Carlo 
procedure, disregarding the terms of order 0{5V/V). 
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for macroscopically planar interfaces. The calculation of the Monte Carlo normalization 
factors does not change the typical scaling of the algorithm, as it consists in calculating the 
histogram for an additional phase of randomly distributed points (which effectively behaves 
as an ideal gas). The better accuracy at larger distances, however, demonstrates that the 
use of the Monte Carlo normalization is much more efficient than the standard approach, as 
it requires much smaller systems to be able to extract the same information (e.g., to resolve 
the fourth peak in Fig. 5, an additional slab of about 2-3 nm would have been needed). 
In this sense, the Monte Carlo normalization procedure can be even beneficial in terms of 
performance. 

V. EXAMPLES OF NON-PLANAR INTERFACES 
A. DPC micelle 

Dodecylphosphocholine (DPC) is a neutral, amphiphilic molecule with a single fatty tail 
that can form micelles in solution: these play a relevant role in biochemistry, especially 
for NMR spectroscopy investigations aiming at understanding the structure of proteins or 
peptides bound to an environment that is similar to the biological membrane^°"^^. The 
molecular structure of DPC is shown in Fig. 6. We have simulated for 500 ps a micelle of 65 
DPC and 6305 water molecules using the force field and configurations from Tieleman and 
collcagucs^^, and have calculated the intrinsic mass density profiles of both phases (DPC 
and water) using GITIM and the Monte Carlo normalization procedure, with a probe sphere 
radius Rp — 0.25. The result of the interfacial atoms identification on the DPC micelle for 
a single frame is shown in Fig. 6, where water molecules have been removed for the sake of 
clarity, and interfacial atoms are highlighted as usual with a halo. The intrinsic mass density 
profile, calculated relative to the DPC surface, is reported in Fig. 7, with the DPC mass 
density profile shown on the left, and the water profile on the right. 

As usual, the delta-like contribution at r = identifies the contribution from interfacial 
DPC atoms. In addition, we have calculated, for the first time, the intrinsic profiles of the 
orientational order parameters Si and 5*2 of the water molecules around the DPC micelle. 
The two parameters are defined as 5"! = (cos(6'i)) and 5*2 (3cos^(6'2) — 1) /2, where 9i and 
02 are the angles between the water molecule position vector (with respect to the micelle 
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center), and the water molecule symmetry axis and molecular plane normal, respectively. 
The orientation is taken so that Oi < 7r/2 when the hydrogen atoms are farther from the 
micelle than the corresponding oxygen. The complete picture of the orientation of water 
molecules would be delivered by the calculation of the probability distribution p{6i, 62)^^'^^, 
but here we limit our analysis to the two separate order parameters and their intrinsic 
profiles. Note that, since these quantities are computed per particle, there is no need to 
apply any volume normalization. The polarization of water molecules, which is proportional 
to Si, appears to be different from zero only very close to the micellar surface. In particular. 
Si has a correlation with the main peak of the intrinsic density profile in the proximity of 
0.4 nm. Water molecules located closer to the interface show a first change in the sign of the 
polarization and a subsequent one when crossing the interface. Farther than 0.25 nm inside 
the micelle, not enough water molecules are found to generate any meaningful statistics. Also 
the S2 order parameter is practically zero beyond 0.6 nm, and again a correlation is seen 
with the main peak of the intrinsic density profile, and the maximum in the orientational 
preference is found just next to the interface, where Si ~ 0, showing that water molecules 
are preferentially laying parallel to the interfacial surface. 

B. Soot 

One of the main byproducts of hydrocarbon flames, soot is thought to have a relevant 
impact on atmospheric chemistry and global surface warming^^'^^. Electron, UV, and atomic 
force microscopy have revealed the size and structure of soot particles from different sources 
at different scales^*^"^^. In particular, soot emitted by aircraft is found to be made of several, 
quasi-spherical, concentric graphitic layers of size in the range from 5 to 50 nm^''. We have 
used four model structures (S]^,S2, S4 and S^^ from Ref. 57) to demonstrate the ability of 
GITIM to identify surface atoms in complex geometries. In Fig. 8, the S\ model is represented 
in section as a triangulated surface (right), showing the four concentric layers, and in whole 
(left) showing the surface atoms as detected by GiTiM using Rp — 0.25 nm. The histograms 
of the total number of atoms and of the surface ones, as a function of the distance from the 
center of the soot particles, are shown in Fig. 9 for the four different models, where it is seen 
how particles of the size of a water molecule have mostly access only to the inner and outer 
parts of the innermost and outermost shell, respectively, and cover them almost completely. 
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This finding is in a clear accordance with the results of the void analysis and adsorption 
isotherm calculations presented in Ref. 57. 

C. Secondary cholic acid micelle 

Bile acids, such as cholic acid are biological amphiphiles built up by a steroid skeleton 
and side groups attached to it. The organization of these side groups is such that hydrophilic 
and hydrophobic groups are located at the two opposite sides of the steroid ring. Thus, bile 
acids have a hydrophilic and a hydrophobic face (often referred to as the a and ^ side, 
respectively) rather than a polar head and an apolar tail, as in the case of other surfactants 
like, for example, DPC. The unusual molecular shape leads to peculiar aggregation behavior 
of bile acids. At relatively low concentrations they form regular micelles with an aggregation 
number of 2-10, while above a second critical micellar concentration these primary micelles 
form larger secondary aggregates by establishing hydrogen bonds between the hydrophilic 
surface groups of the primary micelles^^'^^. These secondary micelles are of rather irregular 
shape,^^'^^ which makes them an excellent test system for our purposes. 

Here we analyze the surface of a secondary cholic acid micelle composed of 35 molecules, 
extracted from a previous simulation work^® and simulated for the present purposes for 500 
ps in aqueous environment. An instantaneous snapshot of the micelle is shown in Fig. 10 
(water molecules are omitted for clarity) together with a schematic structure of the cholic 
acid molecule. We calculated the density profile of water as well as of cholic acid relative to 
the intrinsic surface of the micelle by the GITIM method. The resulting profiles are shown 
in Fig. 11. The micelle has a characteristic elongated shape, which exposes a large part of 
its components to the solvent, so that roughly 80% of the micelle atoms are identified as 
surface ones. The small volume to surface ratio of the micelle is at the origin of the rather 
noisy intrinsic density profile for the micelle itself. The profile, in addition to the delta-like 
contribution at the surface, presents another very sharp peak located at a distance of about 
0.18 nm inside the surface, due to the rather rigid structure of the bile molecule. The water 
intrinsic density profile, on the contrary, shows a marked peak at 0.25 nm, absent in the 
DPC micelle case, due to the presence of hydrogen bonds between water molecules and the 
hydroxyl groups of cholic acid. 
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VI. CONCLUSIONS 



In this paper we presented a new algorithm that combines the advantageous features of 
both the ITIM method^ and the a-shapes algorithm^^'^^ to be used in determining the intrinsic 
surface in molecular simulations. Thus, unlike the original variant, this new, generalized 
version of ITIM, dubbed GITIM, is able to treat interfaces of arbitrary shapes and, at the 
same time, to take into account the excluded volumes of the atoms in the system. It should 
be emphasized that the GiTiM algorithm is not only able to find the external surface of the 
phase of interest, but it also detects the surface of possible internal voids inside the phase. 
The method, based on inflating probe spheres up to a certain radius in points inside the phase 
turned out to provide practically identical results with the original ITIM analysis for planar 
interfaces. Further, its applicability to non-planar interfaces was shown for three previously 
simulated systems, i.e., a quasi-spherical micelle of DPC^^, molecular models of soot^^, and 
a secondary micellar aggregate of irregular shape built up by cholic acid moleculcs^^. 

Another important result of this paper concerns the correct way of calculating density 
proflles relative to intrinsic interfaces, irrespective of whether they are macroscopically planar 
or not. Thus, here we proposed a Monte Carlo-based integration algorithm to estimate the 
volume elements in which points of the proflle are calculated, in order to normalize them 
correctly. The issue of normalization with the volume elements in macroscopically flat fluid 
interfaces originates from the fact that these interfaces are rough on the molecular scale, 
namely, at the length scale of the calculated profiles. We clearly demonstrated that using 
this new normalization the artificial smearing of the intrinsic density profiles far from the 
intrinsic interface can be avoided. 

Two computer programs that implement, respectively, an optimized version of ITIM and 
the new GITIM algorithm, as well as the calculation of intrinsic density and order parameters 
profiles, are made available free of charge at http : / /www . gitim . eu/. The programs are com- 
patible with the trajectory and topology file formats of the Gromacs molecular simulation 
package^^. 
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VII. APPENDIX 

Here, following Ref. 69 we derive the expressions for the radius R and position r = {x, y, z) 
of the center of the sphere which is touching four other ones, having given radii and center 

positions and r-j = {xi,yi,Zi) {i = 1,2,3 or 4), respectively. The conditions of touching 
can be expressed with the following nonlinear system of four equations: 



\r-ri\^ ^{R + Rif. (3) 

By subtracting one of them from the other three (without loss of generality we subtract the 
one with i = 1), the quadratic term, r^, will be eliminated and the system Eq.(3) would 
become linear with respect to r: 

Mr = s - Rd, (4) 

where the matrix M and the vectors d and s are defined as 



ri - r2 
ri - ra 

ri - r4 



d = 



Ri — R2 

Ri — R3 

Ri - Ri 



\ 



(5) 



and 



s = 



rf-rl-Rj + Rl 



\ 



(6) 



ri-vl-Rl + R,! 
rl-rl-Rl + Rl 

Equation (4) has a unique solution if matrix M is non-singular (the singularity of M corre- 
sponds to the case when all 4 spheres are co-planar, which means that the unknown sphere 
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either does not exist, or is not unique): 

r = M-^s - RM-^d = tq - Ru, (7) 

where M~^s = tq and u = M^^d. Once Eq.(7) is substituted into the first of the constraints 
Eq.(3), it leads to the quadratic algebraic equation with respect to R: 

(1 - |u|2) R^ + 2{Ri-u-v)R+{Rl - IvH = 0, (8) 

where v = ri — Tq. The solution of Eq. (8) can be found in the following form: 

-(i^i-u-v)±|i^iu + v| 

^± = T^T^Tp ■ 

If I up is not equal to unity (which corresponds to the case when the 4 spheres are tangential 
to one plane), then Eq.(9) provides two different solutions, and the positive one provides 
the radius R of the touching sphere as a function of the centre position r. Eventually, 
the positions of their centres can be obtained by inserting R into Eq.(7). In the present 
manuscript in case of two possible solutions we choose the sphere with minimal radius. 
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FIG. 1: Left: example of the a-shapes algorithm on a set of points on the plane. The lines connecting 
the atoms represent the Delaunay triangulation (the triangles are labeled by numbers from 1 to 
12). Solid lines mark triangles belonging to the a-complex, and dashed lines those which are not. 
The light-shaded atoms are those belonging to the a-shape, the border of the a-complex. Two 
atoms (in triangle 1) are outside the a-shape, and one (shared by triangles 9-12) is inside the 
a-shape. Right: schematic representation of the itim algorithm, applied to a single water molecule: 
the probe spheres (circles) are moved down the test lines (dashed lines) until they touch an atom. 
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FIG. 2: Simulation snapshot of a H2O/CCI4 system. The oxygen atoms at the interface between the 
H2O phase (inner) and CCI4 phase (outer) as recognized by the gitim algorithm are represented 
with an additional halo. Unconnected points belong to molecules which cross periodic boundary 
conditions. 
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FIG. 3: Average number of surface atoms identified by itim (squares) and gitim (circles) as a 
function of the probe sphere radius. 
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Top view Side view 



FIG. 4: Water surface oxygen atoms in the H2O/CCI4 system in one simulation snapshot as rec- 
ognized by GITIM exclusively (small spheres), itim exclusively (large spheres) or by both methods 
(sphere with halo). 
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FIG. 5: Intrinsic density profiles of water (curves on tlie left) and carbon tetrachloride (curves on 
the right) with respect to the water surface as computed with itim (solid curves) or with GiTiM 
(dashed curves). The thin curves are computed using itim and the Monte Carlo normalization 
procedure described in Sec. IV. 
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FIG. 6: Right: schematic structure of a DPC molecule. Left: snapshot of a DPC micelle in water. 
Only the DPC constituents are shown for the sake of clarity. Atoms with a halo are those recognized 
by GITIM as surface ones. 
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FIG. 7: Upper panel: intrinsic density profiles of water (right) and of DPC (left). Lower panel: 
intrinsic profile of the orientational order parameters (Pi, solid line, P2, dashed line). The vertical 
dashed lines marks the position of the interface. 
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FIG. 8: The S\ soot model^^ represented in section (right, triangulated surface) and in whole 
(left, wireframe) with the atoms identified by gitim as surface ones highlighted using thicker, red 
elements. Besides surface atoms, also chemical bonds between surface atoms are highlighted, as 
well as five, six and seven membered rings (filled surfaces). 
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FIG. 9: Histograms of the atoms in the four soot models taken from Ref. 57. Each panel refers to 
a different structm'e (depicted with wireframe), and presents the distribution of all atoms (filled, 
darker area) and of surface atoms identified by gitim (filled, lighter area), as a function of the 
distance from the center. 



29 




30 




-0.4 -0.2 0.2 0.4 0.6 0.8 
Position relative to the interface / nm 

FIG. 11: Density profile of water (right) and cholic acid (left) in the secondary micelle. 
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